Bogota distance mapping EN

Introduction

This Markdown does the following:

• Mapping the distance.
• Mapping GHG.
• Mapping air pollution

Preliminary work

Defining the working environment

For the moment the data is stored locally. We will further put a link to download the data.

setwd("C:/Users/hugot/Documents/Rennes 2/Calcul distance SIG/Reprise perso/RMD Bogota")

Loading packages

library(sfnetworks)
library(dodgr)
library(sf)
library(sp)
library(spdep)
library(geoR)
library(spatstat)
library(tidygraph)
library(data.table)
library(dplyr)
library(purrr)
library(TSP)
library(openxlsx)
library(raster)
library(stars)
library(terra)
library(ggplot2)
library(ggspatial)
library(knitr)
library(tidyr)
library(cartography)

Loading the data

viajes_dist <- read.xlsx("Tables produites/19-12/ViajesEODH_dist_corr.xlsx")    # Trip table
UTAM <- st_read(dsn = "Data", layer = "EMU2019_area") %>% st_transform(4326)    # UTAM boundaries

## Reading layer `EMU2019_area' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 132 features and 57 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 569984 ymin: 492795.1 xmax: 626083.9 ymax: 557354.6
## Projected CRS: WGS 84 / UTM zone 18N

Manzana <- st_read(dsn = "Data", layer = "Manzana_Hogares")                     # Household manzana

## Reading layer `Manzana_Hogares' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 21828 features and 56 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 569999.6 ymin: 493719.5 xmax: 625763.1 ymax: 557019.5
## Projected CRS: WGS 84 / UTM zone 18N

EF <- read.xlsx("Data/facteurs d'émissions.xlsx")                               # emission factors
limites_bogota <- st_read(dsn = "Data", layer = "Localidad_Municipio_2017") %>% # Municipalities boundaries
  st_transform(4326) 

## Reading layer `Localidad_Municipio_2017' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 564193.1 ymin: 483995.6 xmax: 634113.2 ymax: 570281.8
## Projected CRS: WGS 1984 UTM, Zone 18 North, Meter

Transmi <- st_read(dsn = "Data", layer = "Transmilenio")                        # Transmilenio network

## Reading layer `Transmilenio' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 13 features and 7 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 588057.1 ymin: 503694.6 xmax: 606100.3 ymax: 527398.5
## Projected CRS: WGS 84 / UTM zone 18N

Defining the names to be displayed

centroides_limites_bogota <- st_centroid(limites_bogota) %>%
  mutate(lon = st_coordinates(.)[,1],
         lat = st_coordinates(.)[,2])

centroides_limites_bogota$MUNI <- centroides_limites_bogota$NOMBRE
centroides_limites_bogota$MUNI[centroides_limites_bogota$NOMBR_MUNI == "BOGOTA"] <- ""

centroides_limites_bogota$LOCA <- centroides_limites_bogota$NOMBRE
centroides_limites_bogota$LOCA[centroides_limites_bogota$NOMBR_MUNI != "BOGOTA"] <- ""

centroides_limites_bogota$DISPLAY <- centroides_limites_bogota$NOMBRE
centroides_limites_bogota$DISPLAY[!((centroides_limites_bogota$NOMBR_MUNI == "BOGOTA"  & centroides_limites_bogota$NOMBRE %in% c("BOSA", "SANTA FE", "CIUDAD BOLIVAR", "TEUSAQUILLO", "FONTIBON", "ENGATIVA", "SUBA", "USAQUEN", "USME", "SAN CRISTOBAL", "CHAPINERO", "KENNEDY")) | (centroides_limites_bogota$NOMBR_MUNI != "BOGOTA"))] <- ""

centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "SUBA"] <- 4.719
centroides_limites_bogota$lon[centroides_limites_bogota$NOMBRE == "CHAPINERO"] <- -74.027
centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "USAQUEN"] <- 4.71
centroides_limites_bogota$DISPLAY[centroides_limites_bogota$NOMBRE == "CIUDAD BOLIVAR"] <- "CIUDAD \nBOLIVAR"

Defining the map zoom level and extent

The same zoom level is chosen for Bogotá and Lima for comparison purpose

zoom_level <- 9.5
lon_span <- 360 / 2^zoom_level
lat_span <- 360 / 2^zoom_level

zoom_to_Bogota <- c(-74.12, 4.72)  # Bogota
lon_bounds_Bogota <- c(zoom_to_Bogota[1] - lon_span / 2, zoom_to_Bogota[1] + lon_span / 2)
lat_bounds_Bogota <- c(zoom_to_Bogota[2] - lat_span / 2, zoom_to_Bogota[2] + lat_span / 2)

Distance mapping

Regular mapping

UTAM$Bogota <- 1
contour <- UTAM %>% 
  group_by(Bogota) %>% 
  summarise(geometry = st_union(geometry))

Dist_UTAM <- viajes_dist[,c(1,9,18)] %>%
  left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
  group_by(UTAM)%>%
  summarize(dist = sum(f_exp*dist_real/1000))
  
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")

bks = c(100000,200000,500000,1000000,2500000)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = dist))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#fafaaa", high = "#bb2b30", labels = c("100000", "200000", "500000", "1000000", "2500000"))+
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Total daily PKT per \nUTAM of residence", fill = "Total passengers-km per day", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Dist_UTAM <- viajes_dist[,c(1,6,9,18)] %>%
  filter(p17_Id_motivo_viaje == 1) %>%
  left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
  group_by(UTAM)%>%
  summarize(dist_work = sum(f_exp*dist_real/1000))
  
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")

bks = c(100000,200000,500000,1000000)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = dist_work))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#fafaaa", high = "#bb2b30", labels = c("100000", "200000", "500000", "1000000"))+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Commuting to work", fill = "Passengers-km per day", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Dist_UTAM <- viajes_dist[,c(1,6,9,18)] %>%
  filter(p17_Id_motivo_viaje == 3) %>%
  left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
  group_by(UTAM)%>%
  summarize(dist_school = sum(f_exp*dist_real/1000))
  
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")

bks = c(10000,20000,50000,100000,250000)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = dist_school))+
    geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#fafaaa", high = "#bb2b30") +
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Going to school \nor university", fill = "Passengers-km per day", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Dist_UTAM <- viajes_dist[,c(1,6,9,18)] %>%
  filter(p17_Id_motivo_viaje %in% c(12,16)) %>%
  left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
  group_by(UTAM)%>%
  summarize(dist_leisure = sum(f_exp*dist_real/1000))
  
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")

bks = c(5000,10000,20000,30000,60000)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = dist_leisure))+  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#fafaaa", high = "#bb2b30") +
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Leisure or physical activity", fill = "Passengers-km per day", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Dist_UTAM <- viajes_dist[,c(1,9,10,18)] %>%
  filter(modo_principal == "A pie") %>%
  left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
  group_by(UTAM)%>%
  summarize(dist_walk = sum(f_exp*dist_real/1000))
  
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")

#bks = c(0.5,0.8,1.2,1.5,2)
bks_walk <- round(as.numeric(quantile(UTAM$dist_walk/UTAM$NUMPERSTOT, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- bks_walk

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = dist_walk/NUMPERSTOT))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#f8f8f8", high = "#103D82") +
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Walking", fill = "Average daily distance per capita (km)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Dist_UTAM <- viajes_dist[,c(1,9,10,18)] %>%
  filter(modo_principal == "Auto") %>%
  left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
  group_by(UTAM)%>%
  summarize(dist_auto = sum(f_exp*dist_real/1000))
  
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")

UTAM$dist_auto[UTAM$UTAM %in% c("UTAM3", "UTAM2")] <- NA

#bks = c(1,2,3,4,5)
bks_car <- round(as.numeric(quantile(UTAM$dist_auto/UTAM$NUMPERSTOT, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- bks_car

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = dist_auto/NUMPERSTOT))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#f8f8f8", high = "#103D82", na.value = "grey80") +
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Private car", fill = "Average daily distance per capita (km)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Dist_UTAM <- viajes_dist[,c(1,9,10,18)] %>%
  filter(modo_principal %in% c("Transporte informal", "Bicitaxi")) %>%
  left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
  group_by(UTAM)%>%
  summarize(dist_informal = sum(f_exp*dist_real/1000))
  
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")

bks = c(0.2,0.5,1,2,3)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = dist_informal/NUMPERSTOT))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#f8f8f8", high = "#103D82", na.value = "grey80") +
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Informal transport", fill = "Average daily distance per capita (km)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Dist_UTAM <- viajes_dist[,c(1,9,10,18)] %>%
  filter(modo_principal %in% c("TransMilenio", "Alimentador", "Cable")) %>%
  left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
  group_by(UTAM)%>%
  summarize(dist_BRT = sum(f_exp*dist_real/1000))
  
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")

UTAM$dist_BRT[UTAM$UTAM == "UTAM2"] <- NA

bks = c(1,2,4,6,8)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = dist_BRT/NUMPERSTOT))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_sf(data = Transmi, aes(color = "BRT (Transmilenio)"))+
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#f8f8f8", high = "#103D82",na.value = "grey80") +
  scale_color_manual(values = "red")+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Bus Rapid Transit (BRT)", fill = "Average daily distance per capita (km)", linetype = "", color = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Dist_UTAM <- viajes_dist[,c(1,9,10,18)] %>%
  filter(modo_principal %in% c("TransMilenio", "Alimentador", "Cable", "Intermunicipal", "SITP Provisional", "SITP Zonal")) %>%
  left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
  group_by(UTAM)%>%
  summarize(dist_transit = sum(f_exp*dist_real/1000))
  
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")

UTAM$dist_transit[UTAM$UTAM == "UTAM2"] <- NA

#bks = c(2,4,6,8,10)
bks_transit <- round(as.numeric(quantile(UTAM$dist_transit/UTAM$NUMPERSTOT, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- bks_transit

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = dist_transit/NUMPERSTOT))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#f8f8f8", high = "#103D82",na.value = "grey80")+
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Public transport", fill = "Average daily distance per capita (km)", linetype = "", color = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Spatial smoothing

Preliminary work

UTAM <- UTAM %>% st_transform(32618)
UTAMCentroids <- st_centroid(UTAM)
    
listPPV <- knearneigh(UTAMCentroids, k = 1, longlat = TRUE) # Finding each UTAM nearest neighbor
PPV <- knn2nb(listPPV, row.names = UTAM$UTAM) # Coercing knn objects into nb objects
distPPV <- nbdists(PPV, UTAMCentroids) # Computing the distance between nearest neighbors
print(as.data.frame(t(as.matrix(summary(unlist(distPPV))))))

##       Min. 1st Qu.   Median     Mean  3rd Qu.     Max.
## 1 553.5302 1365.71 1587.616 2181.223 1854.208 10899.47

hist(unlist(distPPV), breaks = 20,
    main = "Distance to the closest neighbor",
    col = "black", border = "white", xlab = "Distance", ylab = "Frequence")

nb <- poly2nb(pl = UTAM, snap = 50, queen = TRUE)
summary(nb)

## Neighbour list object:
## Number of regions: 132 
## Number of nonzero links: 632 
## Percentage nonzero weights: 3.627181 
## Average number of links: 4.787879 
## 16 regions with no links:
## 19 112 113 119 120 121 122 123 124 125 126 128 129 130 131 132
## 18 disjoint connected subgraphs
## Link number distribution:
## 
##  0  1  2  3  4  5  6  7  8  9 10 
## 16  3  4 12 19 19 25 16 11  6  1 
## 3 least connected regions:
## 39 109 111 with 1 link
## 1 most connected region:
## 100 with 10 links

UTAM <- UTAM[!(UTAM$fid %in% c(19, 112, 113, 119, 120, 121, 122, 123, 124, 125, 126, 128, 129, 130, 131, 132)),]
nb <- poly2nb(pl = UTAM, snap = 50, queen = TRUE)
#moran.test(UTAM$dist, listw = nb2listw(nb))

#UTAM$dist_leisure[is.na(UTAM$dist_leisure)] <- 0

moran <- c(0*(1:4))
for(i in 1:4){
m <- moran.test(data.frame(UTAM[,58+i])[,1], listw = nb2listw(nb))
moran[i] <- round(as.numeric(m$estimate[1]), digits = 2)
}

# Displaying the Moran's I
df <- data.frame(cbind(as.numeric(moran),motivo = c("Total", "Trabajo", "Estudio", "Esparcimiento")))
ggplot(df) +  
    theme_classic() +
    geom_point(aes(x = motivo, y=moran), shape = 21, fill = "white", color = "black", size = 2, stroke = 0.5)  +
    labs(x="Trip purpose", y="Moran") +
    theme(legend.position = "left", 
      panel.border = element_rect(colour = "black", fill=NA, linewidth =0.5),
      )

UTAMCentroids <- st_centroid(UTAM)
UTAMCentroids <- as(UTAMCentroids, "Spatial")
UTAMCentroids.geodata <- as.geodata(UTAMCentroids, data.col = "dist")

vario.ex<- variog(UTAMCentroids.geodata, bin.cloud=TRUE, option = "bin")

## variog: computing omnidirectional variogram

plot(vario.ex, main = "Semivariogram of the computed daily distance", cex.main = 1)
lines(vario.ex, type ="l", lty = 2, col="red", add = TRUE)

Smoothed maps

UTAM$Bogota <- 1
contour <- UTAM %>% 
  group_by(Bogota) %>% 
  summarise(geometry = st_union(geometry))

#Creating a color palette
colors <- c("#f9f5a7", "#f7eca0", "#f3d893", "#e9b179", "#d9805a","#c64c3e")
bks = c(0,100000,200000,500000,1000000,2500000,5000000)

# Defining the edge of Bogotá as the extent of the map
Emprise = as.owin(st_buffer(UTAM,0))

# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm) 
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist)

# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))

# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object

# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)

# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages

# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326) 

# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))

ggplot()+
  theme_bw()+
  geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_manual(values = colors, labels = c("< 100000","100000 - 200000","200000 - 500000","500000 - 1000000","1000000 - 2500000","> 2500000"))+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Total daily PKT per \nUTAM of residence", fill = "Total passengers-km per day", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

#Creating a color palette
colors <- c("#f9f5a7", "#f7eca0", "#f3d893", "#e9b179", "#d9805a","#c64c3e")
bks = c(0,20000,50000,100000,200000,500000,2000000)

# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm) 
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_work)

# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))

# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object

# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)

# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages

# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)

# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))

ggplot()+
  theme_bw()+
  geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_manual(values = colors, labels = c("< 20,000","20,000 - 50,000","50,000 - 100,000","100,000 - 200,000","200,000 - 500,000","> 500,000"))+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Commuting to work", fill = "Passengers-km per day", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

#Creating a color palette
colors <- c("#f9f5a7", "#f7eca0", "#f3d893", "#e9b179", "#d9805a","#c64c3e")
bks = c(0,10000,20000,50000,100000,250000,500000)

# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm) 
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_school)

# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))

# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object

# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)

# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages

# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)

# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))

ggplot()+
  theme_bw()+
  geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_manual(values = colors, labels = c("< 10000","10000 - 20000","20000 - 50000","50000 - 100000","100000 - 250000","> 250000"))+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Going to school or university", fill = "Passengers-km per day", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

#Creating a color palette
colors <- c("#f9f5a7", "#f7eca0", "#f3d893", "#e9b179", "#d9805a","#c64c3e")
bks = c(0,5000,10000,20000,30000,60000,100000)

# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm) 
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_leisure)

# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))

# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object

# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)

# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages

# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)

# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))

ggplot()+
  theme_bw()+
  geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_manual(values = colors, labels = c("< 5000","5000 - 10000","10000 - 20000","20000 - 30000","30000 - 60000","> 60000"))+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Leisure or physical activity", fill = "Passengers-km per day", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

UTAM$dist_walk[is.na(UTAM$dist_walk)] <- 0

#Creating a color palette
colors <- c("#efeff3", "#dddee9", "#caccde", "#aeb1ce", "#6c77a9","#284689")
#bks = c(0,0.5,0.8,1.2,1.5,2,4)
bks <- c(bks_walk,100)

# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm) 
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_walk/UTAM$NUMPERSTOT)

# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))

# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object

# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)

# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages

# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)

# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))

ggplot()+
  theme_bw()+
  geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_manual(values = colors, labels = c("< 0.63","0.63 - 0.75","0.75 - 1.03","1.03 - 1.32","1.32 - 1.58","> 1.58"))+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Walking", fill = "Average daily distance per capita (km)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text.align = 1)+   
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+   
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

UTAM$dist_auto[is.na(UTAM$dist_auto)] <- 0

#Creating a color palette
#pal <- colorRampPalette(colors = c("#fafaaa", "#bb2b30"), interpolate = "linear")
#colors <- pal(6) 
colors <- c("#efeff3", "#dddee9", "#caccde", "#aeb1ce", "#6c77a9","#284689")
#bks = c(0,1,2,3,4,5,10)
bks <- c(bks_car,100)

# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm) 
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_auto/UTAM$NUMPERSTOT)

# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))

# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object

# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)

# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages

# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)

# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))

ggplot()+
  theme_bw()+
  geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_manual(values = colors, labels = c("< 0.56","0.56 - 1.01","1.01 - 1.86","1.86 - 3.02","3.02 - 4.65","> 4.65"))+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Private car", fill = "Average daily distance per capita (km)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

UTAM$dist_informal[is.na(UTAM$dist_informal)] <- 0

#Creating a color palette
#pal <- colorRampPalette(colors = c("#fafaaa", "#bb2b30"), interpolate = "linear")
#colors <- pal(6) 
colors <- c("#efeff3", "#dddee9", "#caccde", "#aeb1ce", "#6c77a9","#284689")
bks = c(0,0.2,0.5,1,2,3,10)

# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm) 
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_informal/UTAM$NUMPERSTOT)

# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))


# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object

# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)

# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages

# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)

# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))

ggplot()+
  theme_bw()+
  geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_manual(values = colors, labels = c("< 0.2","0.2 - 0.5","0.5 - 1","1 - 2","2 - 3","> 3"))+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Informal transport", fill = "Average daily distance per capita (km)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

UTAM$dist_BRT[is.na(UTAM$dist_BRT)] <- 0

#Creating a color palette
#pal <- colorRampPalette(colors = c("#fafaaa", "#bb2b30"), interpolate = "linear")
#colors <- pal(6) 
colors <- c("#efeff3", "#dddee9", "#caccde", "#aeb1ce", "#6c77a9","#284689")
bks = c(0,1,2,4,6,8,15)

# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm) 
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_BRT/UTAM$NUMPERSTOT)

# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))


# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object

# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)

# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages

# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)

# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))

ggplot()+
  theme_bw()+
  geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_sf(data = Transmi, aes(color = "BRT (Transmilenio)"))+
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_manual(values = colors, labels = c("< 1","1 - 2","2 - 4","4 - 6","6 - 8","> 8"))+
  scale_color_manual(values = "red")+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Bus Rapid Transit (BRT)", fill = "Average daily distance per capita (km)", linetype = "", color = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

UTAM$dist_transit[is.na(UTAM$dist_transit)] <- 0

#Creating a color palette
#pal <- colorRampPalette(colors = c("#fafaaa", "#bb2b30"), interpolate = "linear")
#colors <- pal(6) 
colors <- c("#efeff3", "#dddee9", "#caccde", "#aeb1ce", "#6c77a9","#284689")
#bks = c(0,2,4,6,8,10,20)
bks <- c(bks_transit,100)

# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm) 
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_transit/UTAM$NUMPERSTOT)

# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))


# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object

# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)

# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages

# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)

# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))

ggplot()+
  theme_bw()+
  geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_manual(values = colors, labels = c("< 3.09","3.09 - 3.97","3.97 - 5.87","5.87 - 7.60","7.60 - 9.24", "> 9.24"))+
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "Public transport", fill = "Average daily distance per capita (km)", linetype = "", color = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Greenhouse gas emissions and air pollutants mapping

Computing emissions

viajes_dist <- viajes_dist %>% left_join(EF, by = "modo_principal")

viajes_dist$`CO2-eq` <- viajes_dist$`CO2-eq`*viajes_dist$dist_real/1000
viajes_dist$CO <- viajes_dist$CO*viajes_dist$dist_real/1000
viajes_dist$NOx <- viajes_dist$NOx*viajes_dist$dist_real/1000
viajes_dist$SO2 <- viajes_dist$SO2*viajes_dist$dist_real/1000
viajes_dist$COV <- viajes_dist$COV*viajes_dist$dist_real/1000
viajes_dist$PM.2.5 <- viajes_dist$PM.2.5*viajes_dist$dist_real/1000
viajes_dist$PM.10 <- viajes_dist$PM.10*viajes_dist$dist_real/1000

Statistics

emissions_mode <- viajes_dist[,c(9,10,18,21:27)] %>%
  group_by(modo_principal) %>%
  summarize(`CO2-eq` = sum(f_exp*`CO2-eq`), CO = sum(f_exp*CO), NOx = sum(f_exp*NOx), SO2 = sum(f_exp*SO2), COV = sum (f_exp*COV), PM.2.5 = sum(f_exp*PM.2.5), PM.10 = sum(f_exp*PM.10))

#CO2 emissions are converted in metric tons.
emissions_mode$`CO2-eq` <- emissions_mode$`CO2-eq`/1000000

#CO emissions are converted in metric tons.
emissions_mode$CO <- emissions_mode$CO/1000000

#NOx emissions are converted in metric tons.
emissions_mode$NOx <- emissions_mode$NOx/1000000

#PM 2.5 emissions are converted in kg.
emissions_mode$PM.2.5 <- emissions_mode$PM.2.5/1000

#CO2 emissions
ggplot(emissions_mode, aes(reorder(modo_principal, -`CO2-eq`), `CO2-eq`))+
  theme_bw()+
  geom_col(fill = 'darkorange2')+
  labs(x="Main mode", y = "Greenhouse gas emissions (tCO2-eq/day)")+
  theme(axis.text.x = element_text(angle = 20, vjust = 0.5, hjust=1))

#CO emissions
ggplot(emissions_mode, aes(reorder(modo_principal, -CO), CO))+
  theme_bw()+
  geom_col(fill = 'darkgreen')+
  labs(x="Main mode", y = "Carbon monoxyde emissions (tCO/day)")+
  theme(axis.text.x = element_text(angle = 20, vjust = 0.5, hjust=1))

#NOx emissions
ggplot(emissions_mode, aes(reorder(modo_principal, -NOx), NOx))+
  theme_bw()+
  geom_col(fill = 'royalblue')+
  labs(x="Main mode", y = "Nitrogen oxydes emissions (t/day)")+
  theme(axis.text.x = element_text(angle = 20, vjust = 0.5, hjust=1))

#PM 2.5 emissions
ggplot(emissions_mode, aes(reorder(modo_principal, -PM.2.5), PM.2.5))+
  theme_bw()+
  geom_col(fill = 'darkred')+
  labs(x="Main mode", y = "PM 2.5 Emissions (kg/day)")+
  theme(axis.text.x = element_text(angle = 20, vjust = 0.5, hjust=1))

Mapping emissions

Total emissions

UTAM <- st_read(dsn = "Data", layer = "EMU2019_area") %>% st_transform(4326)   # UTAM boundaries

## Reading layer `EMU2019_area' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 132 features and 57 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 569984 ymin: 492795.1 xmax: 626083.9 ymax: 557354.6
## Projected CRS: WGS 84 / UTM zone 18N

UTAM$Bogota <- 1
contour <- UTAM %>% 
  group_by(Bogota) %>% 
  summarise(geometry = st_union(geometry))

Emissions_UTAM <- viajes_dist[,c(1,9,21)] %>%
  left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
  group_by(UTAM)%>%
  summarize(CO2 = sum(f_exp*`CO2-eq`/1000000))
  
UTAM <- UTAM %>% left_join(Emissions_UTAM, by = "UTAM")

bks = c(20,50,70,100,150)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = CO2))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#fafaaa", high = "darkorange2")+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "GHG Emissions per UTAM of residence", fill = "Total emissions (tCO2-eq/day)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Emissions_UTAM <- viajes_dist[,c(1,9,22)] %>%
  left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
  group_by(UTAM)%>%
  summarize(CO = sum(f_exp*CO/1000000))
  
UTAM <- UTAM %>% left_join(Emissions_UTAM, by = "UTAM")

bks = c(1,2,3,5,10)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = CO))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#fafaaa", high = "darkgreen")+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "CO Emissions per UTAM of residence", fill = "Total emissions (tCO/day)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Emissions_UTAM <- viajes_dist[,c(1,9,23)] %>%
  left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
  group_by(UTAM)%>%
  summarize(NOx = sum(f_exp*NOx/1000000))
  
UTAM <- UTAM %>% left_join(Emissions_UTAM, by = "UTAM")

bks = c(0.2,0.5,1,1.5,2)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = NOx))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#ebf7ff", high = "royalblue")+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "NOx Emissions per UTAM of residence", fill = "Total emissions (tNOx/day)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Emissions_UTAM <- viajes_dist[,c(1,9,26)] %>%
  left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
  group_by(UTAM)%>%
  summarize(PM.2.5 = sum(f_exp*PM.2.5/1000))
  
UTAM <- UTAM %>% left_join(Emissions_UTAM, by = "UTAM")

bks = c(1,2,5,10,18)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = PM.2.5))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#fafaaa", high = "darkred")+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "PM.2.5 Emissions per UTAM of residence", fill = "Total emissions (kg PM 2.5/day)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Emissions per capita

#bks = c(0.5,0.6,0.8,1,1.2)
bks_co2 <- round(as.numeric(quantile(1000*UTAM$CO2/UTAM$NUMPERSTOT, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- bks_co2

UTAM$CO2[UTAM$UTAM %in% c("UTAM3", "UTAM2")] <- NA

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = 1000*CO2/NUMPERSTOT))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#fafaaa", high = "darkorange2", na.value = "grey80")+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "GHG Emissions per capita \nper UTAM of residence", fill = "Average daily emissions per capita (kg CO2-eq)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

bks = c(30,40,50,60,80)

UTAM$CO[UTAM$UTAM %in% c("UTAM3", "UTAM2")] <- NA

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = 1000000*CO/NUMPERSTOT))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#fafaaa", high = "darkgreen", na.value = "grey80")+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "CO Emissions per capita \nper UTAM of residence", fill = "Average daily emissions per capita (g CO)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

#bks = c(1,2,3,4,6)

UTAM$NOx[UTAM$UTAM %in% c("UTAM3", "UTAM2")] <- NA

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = 1000000*NOx/NUMPERSTOT))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(low = "#ebf7ff", high = "royalblue", na.value = "grey80")+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "NOx Emissions per capita \nper UTAM of residence", fill = "Average daily emissions per capita (g NOx)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

bks = c(60,80,90,120,150)

UTAM$PM.2.5[UTAM$UTAM %in% c("UTAM3", "UTAM2")] <- NA

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = 1000000*PM.2.5/NUMPERSTOT))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_steps(breaks = bks, low = "#fafaaa", high = "darkred", na.value = "grey80")+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "PM.2.5 Emissions per capita \nper UTAM of residence", fill = "Average daily emissions per capita (mg PM 2.5)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

Smoothed emissions per capita

nb <- poly2nb(pl = UTAM, snap = 50, queen = TRUE)
summary(nb)

## Neighbour list object:
## Number of regions: 132 
## Number of nonzero links: 610 
## Percentage nonzero weights: 3.500918 
## Average number of links: 4.621212 
## 16 regions with no links:
## 19 112 113 119 120 121 122 123 124 125 126 128 129 130 131 132
## 18 disjoint connected subgraphs
## Link number distribution:
## 
##  0  1  2  3  4  5  6  7  8  9 10 
## 16  3  4 17 19 18 29  8 11  6  1 
## 3 least connected regions:
## 39 109 111 with 1 link
## 1 most connected region:
## 100 with 10 links

UTAM <- UTAM[!(UTAM$fid %in% c(19, 112, 113, 119, 120, 121, 122, 123, 124, 125, 126, 128, 129, 130, 131, 132)),]
nb <- poly2nb(pl = UTAM, snap = 50, queen = TRUE)
#moran.test(UTAM$dist, listw = nb2listw(nb))

contour <- UTAM %>% 
  group_by(Bogota) %>% 
  summarise(geometry = st_union(geometry))

UTAM$CO2[is.na(UTAM$CO2)] <- 0
UTAM$CO[is.na(UTAM$CO)] <- 0
UTAM$NOx[is.na(UTAM$NOx)] <- 0
UTAM$PM.2.5[is.na(UTAM$PM.2.5)] <- 0

moran <- c(0*(1:4))
for(i in 1:4){
m <- moran.test(data.frame(UTAM[,58+i]/UTAM$NUMPERSTOT)[,1], listw = nb2listw(nb))
moran[i] <- round(as.numeric(m$estimate[1]), digits = 2)
}

# Displaying the Moran's I
df <- data.frame(cbind(as.numeric(moran),motivo = c("GEI", "CO", "NOx", "PM 2.5")))
ggplot(df) +  
    theme_classic() +
    geom_point(aes(x = motivo, y=moran), shape = 21, fill = "white", color = "black", size = 2, stroke = 0.5)  +
    labs(x="Emission type", y="Moran") +
    theme(legend.position = "left", 
      panel.border = element_rect(colour = "black", fill=NA, size=0.5),
      )

#UTAM must now be in projected coordinates
UTAM <- UTAM %>% st_transform(32618)

#Creating a color palette
colors <- c("#FAFAAA", "#F7DF88", "#F5C565", "#F2AA43", "#F09021", "#EE7600")
#bks = c(0,0.5,0.6,0.8,1,1.2,1.5)
bks <- c(bks_co2,50)

# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm) 
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = 1000*UTAM$CO2/UTAM$NUMPERSTOT)

# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))

# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object

# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)

# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages

# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)

# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))

ggplot()+
  theme_bw()+
  geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) + 
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_manual(values = colors, labels = c("< 0.48","0.48 - 0.55","0.55 - 0.69","0.69 - 0.87","0.87 - 1.07","> 1.07"))+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "GHG Emissions per capita", fill = "Average daily emissions per capita (kg CO2-eq)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

#Creating a color palette
#pal <- colorRampPalette(colors = c("#fafaaa", "darkgreen"), interpolate = "linear")
#colors <- pal(6) 
colors <- c("#FAFAAA", "#C8DB88", "#95BE65", "#63A043", "#318221", "#006400")
bks = c(0,30,40,50,60,80,100)

# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm) 
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = 1000000*UTAM$CO/UTAM$NUMPERSTOT)

# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))

# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object

# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)

# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages

# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)

# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))

ggplot()+
  theme_bw()+
  geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_manual(values = colors, labels = c("< 30","30 - 40","40 - 50","50 - 60","60 - 80","> 80"))+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "CO Emissions per capita", fill = "Average daily emissions per capita (g CO)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

#Creating a color palette
#pal <- colorRampPalette(colors = c("#ebf7ff", "royalblue"), interpolate = "linear")
#colors <- pal(6) 
colors <- c("#EBF7FF", "#C9DAF9", "#A6BEF3", "#84A1ED", "#6285E7", "#4169E1")
bks = c(0,4,6,8,10,12,20)

# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm) 
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = 1000000*UTAM$NOx/UTAM$NUMPERSTOT)

# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))

# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object

# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)

# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages

# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)

# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))

ggplot()+
  theme_bw()+
  geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_manual(values = colors, labels = c("< 4","4 - 6","6 - 8","8 - 10","10 - 12","> 12"))+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "NOx Emissions per capita", fill = "Average daily emissions per capita (g NOx)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))

#Creating a color palette
#pal <- colorRampPalette(colors = c("#fafaaa", "darkred"), interpolate = "linear")
#colors <- pal(6) 
colors <- c("#FAFAAA", "#E3C888", "#CD9565", "#B76343", "#A13121", "#8B0000")
bks = c(0,60,80,90,120,150,180)

# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm) 
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = 1000000*UTAM$PM.2.5/UTAM$NUMPERSTOT)

# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))

# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object

# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)

# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages

# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)

# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))

ggplot()+
  theme_bw()+
  geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_sf(data = contour, colour = "grey30", fill = NA) +
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
  scale_fill_manual(values = colors, labels = c("< 60","60 - 80","80 - 90","90 - 120","120 - 150","> 150"))+
  
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(tag = "PM 2.5 Emissions per capita", fill = "Average daily emissions per capita (mg PM 2.5)", linetype = "") +
  theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
        legend.title = element_text(size = 8),
        legend.text = element_text(size=8),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))